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Background: Proton-nucleus collisions have been used as a intermediate baseline for the determi¬ 
nation of cold medium effects. They lie between proton-proton collisions in vacuum and nucleus- 
nucleus collisions which are expected to be dominated by hot matter effects. Modihcations of the 
quark densities in nuclei relative to those of the proton are well established although those of the 
gluons in the nucleus are not well understood. We focus on the effect of these on quarkonium pro¬ 
duction in proton-lead collisions at the LHC at a center of mass energy of 5.02 TeV. 

Purpose: We determine whether it is possible for the LHC proton-lead data to be described by 
nuclear modifications of the parton densities, referred to as shadowing, alone. We compare our 
results to the nuclear modification factor and to the forward-backward ratio, both as a function of 
transverse momentum, pt, and rapidity, y. 

Methods: We employ the color evaporation model of quarkonium production at next-to-leading 
order in the total cross section and leading order in the transverse momentum dependence. We 
use the EPS09 NLO modifications as a standard of comparison. We study the effect of the proton 
parton density and the choice of shadowing parameterization on the pr and rapidity dependence 
of the result. We check whether or not the calculations at leading and next-to-leading order give 
consistent shadowing results. 

We also investigate the size of the mass and scale uncertainties relative to the uncertainty on the 
shadowing parameterization. Finally, we check whether the expected cold matter effect in nucleus- 
nucleus collisions can be modeled as the product of proton-nucleus results at forward and backward 
rapidity. 

Results: We find that the rapidity and pr dependence of the nuclear modihcation factor are gen¬ 
erally consistent with the next-to-leading order calculations in the color evaporation model. The 
forward-backward ratio is more difficult to describe with shadowing alone. The leading and next- 
to-leading order calculations are inconsistent for EPS09 while other available parameterizations are 
consistent. The mass and scale uncertainties on quarkonium production are larger than those of the 
nuclear parton densities. 

Conclusions: While shadowing is consistent with the nuclear suppression factors within the un¬ 
certainties, it is not consistent with the measured forward-backward asymmetry, especially as a 
function of transverse momentum. Data from p -\- p collisions at the same energy are needed. 


I. INTRODUCTION 

In this paper we concentrate on comparison to the 2013 LHC p-|-Pb data on quarkonium production at ,/snn = 5-02 
TeV. The inclusive J/ip and T production data, binned as a function of rapidity and px or rapidity alone (T), 

come from the ALICE [H-Q and LHCb 0,11 collaborations. We briefly discuss the detector acceptances for each 
experiment in turn here but provide a more complete description of the data in Sec. IE 

Runs with mass asymmetric beams are different at the LHC than at RHIC because the LHC beams are not 
symmetric in energy. Instead, in the 2013 LHC p-|-Pb run a 4 TeV proton beam interacted with a 4(Zpb/Apb) = 1.58 
TeV/nucleon Pb beam. The nucleon-nucleon center of mass frame does not coincide with the laboratory frame in an 
energy-asymmetric collision system. Instead the center of mass frame is shifted by Ay = 0.465, taken to be in the 
direction of the proton beam, as defined by experiments. In addition, the ALICE and LHCb detector systems are 
not symmetric around the interaction point, with muon spectrometers on only one side of midrapidity. Therefore the 
beams have to be run in two modes, Pb-|-p and p-|-Pb. In the first case the lead beam is defined to move toward forward 
rapidity while in the second, the proton beam does. Because the second configuration is most similar to fixed-target 
operation and corresponds to small parton momentum fractions in the nucleus, both setups are analyzed according 
to the convention that the proton beam moves to positive rapidity. Thus the case of Pb-|-p collisions corresponds to 
larger parton momentum fractions in the nucleus. 

The ALICE J/ip measurement has been presented in both the central and forward/backward regions. Their muon 
spectrometer covers the pseudorapidity range —4 < yiab < —2.5 in the laboratory frame. Due to the rapidity 
shift in the asymmetric energy beams, the backward rapidity range for dimuon coverage is —4.46 < j/cms < —2.96 
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while the forward rapidity range for dimuon coverage is 2.03 < ?/cms < 3.53. The ALICE central detector includes 
dielectron coverage for \y\a,h < 0.81. The quarkonium measurements in this region correspond to the rapidity range 
-1.37 < j/cms < 0.43. 

The LHCb detector covers 2 < yiab < 5 in the laboratory frame. Like ALICE, they define the forward direction as 
the direction of the proton beam so that their backward coverage is — 5 < j/cms < — 2.5 and the forward coverage is 
1.5 < ycms ^ 4. 

Because there is not yet a measured p + p baseline, the nuclear modification factor. 


RpPh{y,PT) 


dapPb{y,PT)/dydpT 

TpPbd(7pp{y,pT)/dydpT 


( 1 ) 


relies on an interpolation of the p + p cross section. The rapidity shift was taken into account to obtain the p + p 
cross section in the rapidity ranges of the p+Pb measurement. The factor TpPb in the denominator of Eq. ([T|) takes 
the centrality of the collision into account and is calculated in a Glauber model [I|. In this paper, we concentrate 
only on the minimum bias results. See Ref. 0 for a discussion of the centrality dependence of the J/ip measurement 
at RHIC. Work is in progress on the centrality dependence at the LHC [^. 

In addition to studying the nuclear modification factors in these rapidity ranges, the forward-backward ratio, 


Rpsiy^Pr) 


dapphjy > 0,PT)/dydpT 

ddpPbiy < 0,pT)/dydpT 


( 2 ) 


has also been presented Si- The unmeasured quantities TpPb and app cancel for a rapidity range symmetric around 
2/cms = 0. Therefore Rps is formed in the rapidity region where the forward and backward acceptances completely 
overlap, 2.96 < |ycms| < 3.53 for ALICE and 2.5 < |?/cms| < 4 for LHCb. Some systematic uncertainties also cancel 
in the ratio. This ratio is perhaps a more faithful representation of cold matter effects on the p-l-Pb cross section. 
However, for theoretical interpretation, Rpph is still desirable because even a wrong model can produce the right ratio. 

In a previous paper Q, predictions were made for pA collisions at y/s^ = 8.8 and 5.5 TeV and ratios were formed 
both to p + p collisions at the same energy and to p + p collisions at the anticipated top energy of ^/s = 14 TeV. 
These calculations were made before the LHC turned on and so did not employ the same energies at which data were 
ultimately taken in LHC Run I: ^ = 2.76, 7 and 8 TeV for p+p collisions and 5.02 TeV for p-|-Pb. These calculations 
assumed that the leading and next-to-leading order treatments of the modifications of the parton densities in nuclei, 
when employed consistently, would be identical, as discussed in more detail later. 

More recently, calculations were made for the nuclear modification factor as a function of y and pt at the energy 
appropriate for the p-|-Pb run Q but not taking the rapidity shift into account for the pt acceptance. In addition, 
the Pt dependent ratio was presented for forward rapidity only. These predictions, along with other, updated, 
CEM calculations with EPS09 NLO nuclear parton densities (nPDEs) were compared to the ALICE and LHCb J/ip 
and T data EMi, [iM3- In those calculations, the incorrect factorization scale, pp-: was passed to the nuclear 
parton densities. (The square of the scale was passed instead of the scale itself, as required by most shadowing 
parameterizations, so that the overall cold matter effect was reduced relative to the true value. This was not the case 
for the LO predictions in Ref. Q.) In this paper, this error is corrected and the pr-dependent ratios calculated in the 
CEM are presented for the first time. 

In the next section. Sec. m we provide a brief summary of the ALICE and LHCb J/^/> and T measurements to 
place our calculations in context. 

We then summarize our calculation of quarkonium production in p -|- p collisions in Sec. m We compare the J14’ 
and T distributions obtained with several different sets of proton parton densities. Comparison of these calculations 
to available p -I- p data can be found in Refs, [ij, [l^. A short summary of cold nuclear matter effects is given in 
Sec.HVl 

The nuclear parton density modifications used in this paper are described in Sec. El The calculations are compared 
to the p-l-Pb data from ALICE and LHCb on J/4 and T production at = 5 TeV in Sec. ED We calculate 

both RpPh and Rfb as functions of rapidity and transverse momentum. The data are first compared to the EPS09 
uncertainty bands. We contrast the leading order (LO) and next-to-leading order (NLO) results for order-by-order 
consistency. Next, the data are compared to all the nuclear parton densities discussed in Sec. El to see if any 
parameterizations are particularly favored. The mass and scale dependence of the results is also shown. Finally, we 
discuss how closely the A + A calculations can be reproduced by a convolution of p -I- A and A -|- p collisions and also 
compare to the RHIC data. 
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II. DESCRIPTION OF THE p+PB QUARKONIUM DATA 

In this section, we describe the quarkonium data for J/il) and T(15') from ALICE and LHCb. Since there has been 
no p + p run at ^/s = 5 TeV to date, the denominator of i?ppb has to be interpolated between available measurements 
at other energies. The p + p cross sections also had to be adjusted to the rapidity ranges of the p+Pb measurement. 

The interpolation methods used by the collaborations depend on the quarkonium state, the observable, and the 
previously available data. However, in all cases, the quarkonium states were assumed to be produced unpolarized. 

To obtain RpPh{y) for the J/V'j the ALICE Collaboration used an energy interpolation between their p + p mea¬ 
surements at 2.76 and 7 TeV to obtain the ^/s dependence of the p + p cross section. They presented their J/%1) p+Pb 
data in both a single rapidity interval as well as in six rapidity bins measured for p + p collisions (of course without 
the rapidity shift in p + p). These first results were for i?ppb at forward and backward rapidity and Rfb as a function 
of pt and y [l| . 

Their energy interpolation was based on three assumed shapes: linear, power law and exponential. The central 
value of the result for each rapidity bin is an average of the three shapes while the uncertainty is the quadrature sum 
of a term related to the uncertainty on the data used for the interpolation and a term related to the spread between 
results with different shapes. 

An additional small systematic uncertainty was obtained by comparing the shapes with those of the leading order 
CEM and the fixed-order next-to-leading logarithm (FONLL) approach for inclusive open heavy flavor production [l|. 
We note that the LO and NLO CEM energy dependence should be similar if the same mass and scale parameters, 
as well as the same proton parton densities are used. (The shapes will not be similar if e.g. CTEQ6M is used at 
NLO and CTEQ61L is used at LO.) Using the FONLL calculation for the total cc cross section may produce a shape 
similar to the CEM but the magnitude, of course, will be quite different. 

In a later paper, the ALICE Collaboration presented results for the midrapidity i?ppb, a bin around —1.37 < 
Veins < —0.43 to add to RpPhiv) as well as the ratios RpPhiPr) at forward, backward and midrapidity [^. While the 
forward-backward ratio as a function of pt was published in Ref. [I|, the separate values of RpPhiPr) were not yet 
available. 

The p + p baseline for the pr-dependent ratios was obtained through interpolation. At midrapidity, data from 
i/s = 0.2, 1.96, 2.76 and 7 TeV were used. The 1.96 TeV results from pp collisions from the Tevatron were considered 
on the same basis as p + p collisions because, at these high energies, production is dominated by the gg process. 
Scaling in xt = rriTjy/s was used to compare the disparate energies. Only exponential, logarithmic and power law 
dependencies were considered in this case because there is no pt dependence in the LO CEM and the FONLL approach 
is for single inclusive distributions, not pairs, so the px slope is not available from these calculations. 

At forward rapidity, the only data available to include in the y/s interpolation are the 2.76 and 7 TeV data from 
ALICE. (They did not employ the LHCb results in their interpolation.) In this region, the dependencies were linear, 
power law and exponential. The results for RpPhiPx) are only shown up to 8 GeV because the p + p data were limited 
to this Pt ra ng e. The non-prompt production from b decays increases with pr, giving a ^ 20% correction at 
Pt^8 GeV 0. 

ALICE has also measured the ijj' inp+Pb collisions, finding significantly more suppression [l5j |. Since this difference 
cannot be due to initial state effects alone, we do not address that result in this work. 

The ALICE Collaboration also measured the inclusive T(IS) and T(2S) rates. The rapidity dependence of RpPh 
was reported in Ref. 0. The T yields are not large so that only one rapidity bin is reported. To obtain the p + p 
baseline for RpPh, they used the LHCb results for T production at i/s = 2.76, 7 and 8 TeV, divided into rapidity bins. 
They employed 21 different shapes, 15 from the LO CEM with different proton PDFs and factorization scale choices; 
3 based on the FONLL h quark distributions; while linear, power-law, and exponential shapes rounded out the set. 
The agreement of all the shapes with the data was generally poor so the fits with the worst x^/dof were discarded for 
the final fits. In addition to the J/V’ uncertainties described above, they also considered small rapidity shifts between 
the ALICE and LHCb rapidity bins. 

They found that the T Rpp^ is quite similar to that of J/tp at forward rapidity while at negative rapidity the T 
RpPh is compatible with unity but lower than that of the J/ip. This is a fairly remarkable result since nuclear effects 
are generally expected to be smaller for the T than for J/ip so that the T RpPh should be closer to unity at low px 
where the difference in mass is the dominant effect. At higher px, where px ^ rn, the results should be similar for 
the two quarkonium states. 

They reported the T(2S)/T(1S) ratios at forward and backward rapidity, 0.26 ± 0.09 ± 0.04 and 0.27 ± 0.08 ± 0.04 
respectively 0 . The result is consistent with the p + p ratio at 7 TeV. This is also consistent with shadowing being 


^ Note that the FONLL approach calculates the single inclusive heavy flavor distributions, not those of the QQ pair as in the CEM. 
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Experiment 

y acceptance 

Nj/ip [Ref.] 

Nr(is) [Ref.] 

Nr( 2 S) [Ref.] 

Vt( 3S) [Ref] 

ALICE 

2.03 < j/cms < 3.53 

(6.69 ± 0.05) X 10‘‘ [1] 

305 ± 34 [3j 

83 ± 23 [3j 

- 


—4.46 < 2/cms < —2.96 

(5.67 ±0.05) X 10"^ [1] 

161 ±21 [3j 

42 ± 14 [3j 

- 


-1.37 < ycms < -0.43 

465 ± 37 [2j 

- 

- 

- 

LHCb 

2.5 < Perns < 4.0 

25280 ± 240 [Tj 

189 ± 16 

41 ±9 

13±7 


-4.0 < Perns < -2.5 

8830 ± 160 \£ 

72 ± 14 

17± 10 

4±8 


TABLE I: The J/tp and T yields from the ALICE and LHCb collaborations in their stated rapidity acceptance. The results 
are integrated over all pr ■ 


the dominant cold matter effect since it affects the excited states the same way as the ground state. 

The LHCb Collaboration has also measured 0 and T [1| production in their muon spectrometer. In addition 
to the inclusive J/ip result, they also separate b ^ J/ip decays to present a non-prompt J/ip result. Their primary 
functional dependence to interpolate between their p+p results at 2.76, 7 and 8 TeV is a power law, a{^/s) = (v^/po)^^ ■ 
They use linear and exponential dependencies to obtain a systematic uncertainty on their interpolation. They do not 
use any production models for the energy interpolation [^. A similar, power-law-based method is used to obtain the 
p + p baseline for T production Q . 

Recently, the ATLAS Collaboration has presented results on the forward-backward J/'0 ratio as a function of y in 
the region 8 < pt < 30 GeV and as a function of rapidity in the rapidity range |?/cms| < 1-94 [I^. Their results for 
Rfb are consistent with unity within the uncertainties of the data. CMS also recently presented the prompt J/ip 
Rfb as a function of pp and as a function of rapidity for high pp [a. Their results are consistent with those of 
ATLAS. 

The LHCb T measurement includes low statistics for the T(3S) states, as well as the IS and 2S states. The ratios 
T(2S)/T(1S) are 0.28±0.14±0.05 in the backward direction and 0.20±0.05±0.01 at forward rapidity, both consistent 
with p-bp measurements While the T(3S)/T(1S) ratios are also consistent with those in p-bp collisions, their 
low statistics reduces their significance. 

Finally, we note that we do not discuss the intriguing CMS T(nS)/T(lS) ratios measured as a function of both the 
number of tracks at midrapidity and the transverse energy at forward rapidity M- This effect, also seen in p -bp 
collisions, is not attributable to initial-state modifications of the parton densities and, as such, is outside the scope of 
this work. 

The quarkonium yields for all the ALICE and LHCb results discussed here are given in Table ID 


III. PRODUCTION IN p -b p COLLISIONS 

Following our previous work [ 1 , we treat quarkonium production within the color evaporation model (CEM). In 
the CEM, heavy flavor and quarkonium production are treated on an equal footing. The CEM has enjoyed considerable 
phenomenological success when applied at NLO in the total cross section and LO in the quarkonium pp distribution 
|l^ [igl - l^ . (See Ref. [l^ for comparison of the y/s = 2.76 and 7 TeV ALICE data with the same CEM calculation 
employed here.) 


A. Color Evaporation Model Calculation 


In the CEM, the quarkonium production cross section is some fraction, Fq, of all QQ pairs below the HH threshold 
where H is the lowest mass heavy-flavor hadron. Thus the CEM cross section is simply the QQ production cross 
section with a cut on the pair mass but without any constraints on the color or spin of the final state. The color 
of the octet QQ state is ‘evaporated’ through an unspecified process which does not change the momentum. The 
additional energy needed to produce heavy-flavored hadrons when the partonic center-of-mass energy, yfi, is less than 
2mH, the HH threshold energy, is nonperturbatively obtained from the color field in the interaction region. Thus 
the quarkonium yield may be only a small fraction of the total QQ cross section below 2mH- At leading order, the 
production cross section of quarkonium state C in a p -b p collision is 


^CEM 


)=FcJ2 


pAmi 


4m2 


ds 


dxidx2 f^{xi,pp) fj{x2,p%) J{s)d^j{s,p%,pp) , 


( 3 ) 
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where ij = qq or gg and o'ij(s) is the ij —?> QQ subprocess cross section. Here J'(s) is a kinematics-dependent 
Jacobian. At LO J/(s) = S(s — xiX 2 s)/s, at NLO and for differential cross sections, the expressions are more complex. 

The fraction Fc must be universal so that, once it is fixed by data, the quarkonium production ratios should be 
constant as a function of ■\/s, y and The actual value of Fq depends on the heavy quark mass, m, the scale 
parameters, the parton densities and the order of the calculation. 

We fit Fj!^ to both the full data set as well as to more limited sets. Our final J/'0 result is based on the total cross 
section data with only p, Be, Li, C, and Si targets respectively. In this way, we avoid uncertainties due to ignoring 
any cold nuclear matter effects which are on the order of a few percent in light targets. We also restricted ourselves 
to the forward cross sections only, rather than include the Buda/dy\y=Q data in the fits. The rapidity distributions 
calculated in the MNR code are subject to fluctuations about the mean, even with high statistics calculations. The 
total cross sections, not subject to these fluctuations, are thus more accurate. See Ref. [I^ for more detail. 

In the case of T production, however, most of the reported cross section values are for Buda/dy\y=Q. The branching 
ratio Bii here is a composite for the three T S states which were not separated at fixed-target energies. For later 
experiments, with sufficient resolution to separate the mass peaks, the individual y = 0 cross sections were multiplied 
by the PDG values of the branching ratios and summed [IJ. The data in the T fits are from fixed-target energies, 
19.4 < y/s < 44 GeV, and collider data from the ISR, SppS and the Tevatron. The p + p data from the SppS and the 
Tevatron are fit with the same coefficient as the lower energy p + p data. At ^/s = 630 GeV, the difference between 
the p + p and p + p bb cross sections is less than 0.5%, too small to affect the fit results. 

We use the same values of the charm quark mass and scale parameters as found in Ref. [l^ to obtain the normal¬ 
ization Fc for the J/V', {m, pp/m, pn/m) = (1.27 ± 0.09 GeV, 2.1^Qgg, 1.6^° In the case of T production, we 
employ (m, pp/in, ppi/m) = (4.65 ± 0.09 GeV, l.lj'jg'^o)- We determine Fq only for the central parameter set 

in each case and scale all the other calculations by the same value of Fc to obtain the extent of the J/'!/’ and T mass 
and scale uncertainty bands, as described in detail in Sec. I VI El 

We find Fjf^ = 0.020393 for the central result with {m, pp /m, pR/m) = (1.27 GeV, 2.1,1.6) employing the GTIO 
parton densities dl. We obtain = 0.0077 for the central y = 0 result for the combined T S states with the 
CTIO parton densities and {m, pp/m, pR/m) = (4.65 GeV, 1.4,1.1). After separating the S states, the inclusive IF 
value is Fy — 0.022 [T3| . 

Our GEM calculations use the NLO QQ code of Mangano et al. (MNR) [2^ with the FI FI mass cut in Eq. ([3|), 
as described in Refs. dilii. Because the NLO QQ code is an exclusive calculation, we take the mass cut on the 
invariant average over kinematic variables of the c and c. Thus, instead of defining pp and pR relative to the quark 
mass, they are defined relative to the transverse mass, pp,r oc rriT = \/ rn? + p^ where pp is that of the QQ pair, 
p\ = 0.5(p|,^ APt-)- 

At LO in the total cross section, the QQ pair pp is zero. Thus, while our calculation is NLO in the total cross 
section, it is LO in the quarkonium px distributions. In the exclusive NLO calculation both the Q and Q variables 
are integrated to obtain the pair distributions, recall pp.r oc mp. 

Results on open heavy flavors indicate that some level of transverse momentum broadening is needed to obtain 
agreement with the low pp data. This is often done by including some intrinsic transverse momentum, kp, smearing to 
the initial-state parton densities. The implementation of intrinsic kp in the MNR code is not handled in the same way 
as calculations of other hard processes due to the nature of the code. In the MNR code, the cancellation of divergences 
is done numerically. Since adding additional numerical Monte-Carlo integrations would slow the simulation of events, 
in addition to requiring multiple runs with the same setup but different intrinsic kp kicks, the kick is added in the 
final, rather than the initial, state. In Eq. ®, the Gaussian function gp(kp), 

HpikT) = exp(-fc|/(/c|,)p) , (4) 

[ 2 ^, multiplies the parton distribution functions for both hadrons, assuming the x and kp dependencies in the initial 
partons completely factorize. If factorization applies, it does not matter whether the kp dependence appears in the 
initial or final state if the kick is not too large. In Ref. [ 2 ^, {k^)p = 1 GeV^, along with the Peterson fragmentation 
function with parameter e = 0.06, was found to best describe fixed-target charm production. We note that currently 
Peterson fragmentation with e = 0.06 is considered too strong. The EONLL fragmentation scheme for open heavy 
flavor is softer [ 2 ^. 

In the code, the QQ system is boosted to rest from its longitudinal center-of-mass frame. Intrinsic transverse 
momenta of the incoming partons, kpi and kp2, are chosen at random with k^^ and k‘^2 distributed according to 
Eq. (IH). A second transverse boost out of the pair rest frame changes the initial transverse momentum of the QQ 
pair, Pp, to pp + kpi J- kp2. The initial kp of the partons could have alternatively been given to the entire final-state 
system, as is essentially done if applied in the initial state, instead of to the QQ pair. There is no difference if the 
calculation is LO but at NLO an additional light parton can also appear in the final state so the correspondence is 
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not exact. In Ref. [2^, the difference between the two implementations is claimed to be small if < 2 GeV^. We 
note that the rapidity distribution, integrated over all px, is unaffected by the intrinsic kx- 

The effect of the intrinsic kx on the shape of the quarkonium px distribution can be expected to decrease as ^/s 
increases because the average px also increases with energy. However, the value of may increase with ^/s so that 
effect remains important at higher energies. We assume the form = 1 + (1/n) ln(-ys/20) GeV^. Using the RHIC 
J/U data, we found that n = 12 gave the best description of the px distribution both at central and forward rapidity 
[ISj. A larger value of n and thus of (fcy) is required for the T px distribution. We set n = 3 for T from comparison 
to the Tevatron results at i/s = 1.8 TeV M- For this study, n is not modified by the nuclear medium. 

We note that most approaches to quarkonium production: the GEM; the color singlet model (GSM); and the 
Nonrelativistic QGD approach (NRQGD), assume the validity of collinear factorization which separates the initial, 
nonperturbative parton densities from the perturbatively-calculable hard scattering that produces the final state. We 
assume collinear factorization to hold for quarkonium production in the GEM. Factorization has been proved for 
quarkonium production in NRQGD at high px but not at low px- The open charm flavor results at low px 
at the LHG fisl agree with calculations employing collinear factorization [l3l| better than calculations employing kx 
factorization [27j| . Collinear factorization should work better for bottom since the factorization scale and the x region 
probed are both larger. 

Since the start up of the LHC, several groups have performed global analyses of the nonperturbative matrix elements 
in the NRQGD approach up to NLO, see Ref. [2^ and references therein for discussion and comparison of the results 
for J/ip production data in e+ + e“, e + p, p + p and p + p collisions. In Ref. [2^, it is clear that the nonperturbative 
matrix elements are quite sensitive to the minimum px employed in the fits. These matrix elements do not appear to 
be universal since choosing different data sets to fit to result in quite different values of the matrix elements. Indeed, 
the agreement with the e’*' + e~ and e + p data is poor unless the minimum px is low, px 3 GeV, but these fits 
cannot reproduce the measured high px quarkonium polarization [2^ . 

In addition, if the fitted matrix elements are used to calculate the J/iIj cross section at y = 0 in p + p collisions 
as a function of y/s, good agreement with neither the shape nor the magnitude of the cross section can be obtained, 
see Ref. [2^. The results overshoot the measured cross sections significantly, sometimes by an order of magnitude 
[^ . This is not a surprising outcome because the the integrated y = 0 cross section is dominated by low px J/'0 
production. In the same paper, the GSM cross sections are also compared to these data. It was found that only the 
LO GSM calculation produces a physical ^/s dependence that agrees relatively well with the data, at NLO some values 
oi fix give an unphysical energy dependence [^. They also point out that the GEM produces the best agreement 
with the dependence over the entire range. 

We also note that, while the GSM and NRQGD make predictions for the px dependence at relatively high px, the 
region of interest for nuclear effects on the parton densities, as addressed here, is at low px and the GEM, through 
the application of kx smearing, is the only production model that addresses the entire px range. Since the average 
kx is an energy-dependent parameter, ideally this should be replaced by a low-p^ resummation. 

Some recent results may provide improvements for NRQGD, at least at collider energies [2^, HO,[Ml. References [s^, 
Hll perform small x resummation in the color-glass condensate (CGG) in the NLO NRQGD approach with the 
nonperturbative matrix elements taken from the high px fits of Chao et al. (3^ . The CGG gluon distribution is 
employed in the CGC-I-NRQCD calculation at low px- This matches well with the high px NLO NRQGD calculation 
in the intermediate px range. Thus the entire px range of the data can be described within the combined approaches, 
both in p + p and p + A collisions for forward rapidities and sufficiently high [13, [HI. The energy dependence 
of da/dy\y=o is also reproduced quite well for y/s > 0.2 TeV athough agreement with the midrapidity RHIC data is 
rather poor [30j| - This approach is inapplicable at lower i/s. Reference [26l| presents a factorized power expansion 
for quarkonium production, including next-to-leading power (NLP) contributions to the perturbative part. This 
formalism requires fragmentation functions for heavy QQ pairs as well as for light partons. With the fragmentation 
functions and the NLP contributions included in NRQGD, it was found that the and components of the 
cross section are dominated by the NLP contributions over allpx, independent of the nonperturbative matrix elements 
[131 • This formalism describes J/ip production well at collider energies for px > 10 GeV. The possible dominance 
of the contribution in the total production rate could explain the apparent unpolarized J/'ij) production. This 
conclusion is consistent with the data-driven approach to polarization in Ref. [33l| . 

We note that the polarization has not yet been calculated in the GEM. While it should be straightforward at LO, 
to go to NLO one would have to start from a polarized QQ pair production code. 


B. Comparison ot p + p Results 

We now turn to a comparison of the J and T px and rapidity distribution in p + p collisions for different proton 
parton densities. Our main results are obtained with the GTIO |34l | parton densities, also used in our evaluations of 
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FIG. 1: (Color online) The gluon distribution function at the factorization scale for J/'ip (a) and T (b) production. The CTIO 
(black), CTEQ5M (blue), CTEQ6M (red) and MSTW (magenta) are compared, all calculated with the same input parameters. 


the charm and bottom cross sections MM where we find good agreement with the present quarkonium data within 
the mass and scale uncertainties. In Fig. [1] we compare the NLO CTIO gluon distributions to those from CTEQ5M 
[ 3 ^, CTEQ6M [ 3 ^ and MSTW [13 . The MSTW central NLO set is a rather recent, frequently used parton density 
shown as an alternative to CTIO. It is available at LO, NLO and NNLO and has a low starting scale of /Tq = 1 GeV^. 
The CTEQ sets and CTIO, at LO and NLO, have a starting scale of fig = X GeV^. The vertical lines at the top of 
the plots indicate the gluon density at forward, mid- and backward rapidity. 

The CTEQ6M and GTEQ5M distributions were chosen because they have been used to extract the modifications 
of the parton densities in the nucleus. The CTEQ6M distributions were used in a global analysis of the nuclear 
modifications of the parton densities to obtain the EPS09 NLO sets [13 • They were also used in previous estimates 
of the quarkonium cross sections at the LHC [13 . The corresponding EPS09 LO sets were obtained based on the 
CTEQ61L densities, as we discuss in more detail later. The CTEQ5M distributions were used to extract the NLO 
nuclear modifications in the FGS sets [dOj. Older distributions used to extract nuclear modifications are CTEQ4L 
[m and GRV LO [13 (EKS98 0,113) and GRV98 LO and NLO [13 (nDS [13)- Since these proton PDFs are now 
outdated, we do not show them here. 

Figure [1] shows the gluon densities in the proton at the scales used to calculate J/ijj (a) and T (b) production. The 
CTIO distributions, the most recent of all those considered, follow the previous CTEQ6 global analysis. The CTIO 
and CTEQ6M gluon distributions are almost identical. The earlier CTEQ5M set is quite different: the lowest x value 
included in the analysis is I0“® and, instead of extrapolating smoothly to lower values of x, the density is frozen such 
that xg{x < 10“®, = xg{x = 10“®, It also tends to be higher than the other gluon PDFs for 10“® < x < 10“^. 

The MSTW distribution is below the others for the J/ijj factorization scale. All the gluon distributions are more 
similar at the T scale even though CTEQ5M is still somewhat higher for I0“® < x < 10“^. 

The CTEQ6M and CTIO gluon distributions are zero at the minimum scale, xg{x,gl) = 0. Thus they undergo 
rapid evolution at low x since they are based on a valence-like initial distribution, a;“(l — x)^. On the other hand, 
the MSTW sets allow a negative gluon density at low x. Therefore the MSTW gluon density exhibits a slower scale 
evolution. It will thus produce the smallest cross sections for J/ip and T production while CTEQ5M will give the 
largest. 

In Fig. we show the shapes of the J/ip pr and y distributions for the proton parton densities presented in Fig. [T] 
The pt distributions are shown in the region of rapidity overlap with the shifted p-l-Pb range at forward rapidity, 
2.96 < y < 3.53. In p-\- p collisions, the pt distributions at forward and backward rapidity are identical because the 
y distributions are symmetric. 

The results are given with the same CEM normalization, Fq, as found for the CTIO fits in Ref. [l3. Because the 
low X gluon distributions differ in shape and magnitude, the values of the cross sections also differ by as much as 
15 — 20%. If fits of Fc were made with the three additional proton parton densities, the values of Fc would clearly 
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FIG. 2: (Color online) The J/ij) px distribution at forward rapidity in p + p collisions (a) and the pT-integrated y distribution 
(b). Results from CTIO (black), CTEQ5M (blue), CTEQ6M (red) and MSTW (magenta) are compared, all calculated with 
the same input parameters and using the same value of Fc as for CTIO. 
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FIG. 3: (Color online) The ratios of J/'4> production as a function of pr at forward rapidity in p + p collisions (a) and y 
integrated over all px (b) relative to CTIO for CTEQ5M (blue), CTEQ6M (red) and MSTW (magenta). 


differ. By employing the same value of Fc for all the sets, we emphasize the differences in magnitude as well as shape. 

The CTEQ5M distribution, higher than the other gluon densities at the J/i/i scale for x < 0.01, results in a 20% 
larger overall cross section. This manifests itself at low px where the integral difference is 32% and, most prominently, 
in the rapidity distribution around y ^ 0 where it is ^ 50% higher. The CTEQ5M gluon distribution, as already 
noted, takes a constant value for x < 10“®. The rapidity distribution for the PDE is consequently narrower than the 
other three shown. The corresponding pp distributions do not reflect the x < 10“® behavior because, in the rapidity 
range illustrated, the low pp J/'ijj's are at higher x values, see the y = 4 line in Fig. [T](a). 

The similarity of the CTIO and CTEQ6M gluon distributions results in very similar pp and y distributions. The 
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FIG. 4: (Color online) The T distribution at forward rapidity in p + p collisions (a) and the pr-integrated y distribution (b). 
Results from CTIO (black), CTEQ5M (blue), CTEQ6M (red) and MSTW (magenta) are compared, all calculated with the 
same input parameters and using the same value of Fc as for CTIO. 


MSTW set gives a 15% lower cross section for the J/ip with most of the difference at lowpp. The rapidity distributions 
are also slightly narrower than CTIO but less than that of the CTEQ5M calculation. 

Figure [3] presents the ratios of the CTEQ5M, CTEQ6M and MSTW calculations to the CTIO results as a function 
of pt and y. Aside from the lowest pr bins, the ratios in Fig. [SKa) are practically independent of pr except for 
CTEQ5M which seems to show slow growth with pr after a strong drop at pp < 5 GeV. For 5 < pp < 30 GeV, 
the ratios to CTIO all agree within 10%. The CTEQ6M result differs from CTIO by ^ 2% for pp > 5 GeV. The 
pp-integrated rapidity distributions in Fig. |3] are all quite different. While the GTEQ6M result is within 5% of that of 
GTIO over all rapidity, the other ratios are narrower. Note that the pp > 0 ratio in Fig. [3Kb) is larger than the ratio 
at Pp > 5 GeV in|3Ka) because the y distribution is dominated by low pp. The sharp drop for CTEQ5M is where 
X ^ 10“^. Since the MSTW gluon distributions do not have the same behavior, they are broader than GTEQ5M but 
the ratio is less than 0.7. 

The corresponding T results are displayed in Figs. |4| and |5j The larger factorization scale makes the results all very 
similar in magnitude for this case with the integrated cross sections differing by less than 10%, even for CTEQ5M. 
The ratios as a function of pp are all relatively flat, especially between GTEQ6M and GTIO. The largest difference 
among the results appears at large rapidities where there is a sudden increase in the ratio for |y| ~ 5, dropping off 
afterward where x ^ 10“^ is reached for T production. 

We have shown that the rapidity distributions are more dependent on the chosen baseline proton parton densities 
than the pp distributions. This arises because of the differing behavior of the low x gluon densities. Due to the larger 
scale, the dependence on the proton parton densities is less important for the higher mass T. 


IV. COLD NUCLEAR MATTER EFFECTS 

To go beyond p + p collisions, the proton parton densities must be replaced by those of the nucleus. We assume 
that if A is a nucleus, the nuclear parton densities, ff^{x 2 , factorize into the nucleon parton density, ff{x 2 ,p^), 
independent of A; and a shadowing ratio, Sp g(A, that parameterizes the modifications of the nucleon parton 

densities in the nucleus. The centrality dependence will be discussed elsewhere [ 3 ] so here we deal with S^{A,X 2 ,p'^) 
alone. (Note that we refer to X 2 here since this is the 

We do not consider other cold matter effects in this paper because we want to determine how far the data can be 
described by the assumption of shadowing alone. However, since lower energy J/ij} production clearly requires some 
effects beyond shadowing, we briefly mention these here. 

At fixed-target energies, the xp dependence clearly shows that shadowing is not the only contribution to the 
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FIG. 5: (Color online) The ratios of T production as a function of pr at forward rapidity in p + p collisions (a) and y integrated 
over all pr (b) relative to CTIO for CTEQ5M (blue), CTEQ6M (red) and MSTW (magenta). 


J/ijj nuclear dependence as a function of Feynman xf = 2pz/^/s for the final-state particle [ 13 , Indeed, the 
characteristic decrease of ai^xp) (in dapA/dxp = [dupp/dxF)A°‘^^’^^) for xp > 0.25 cannot be explained by shadowing 
alone In fact, the data so far suggest approximate scaling with xp, not the target momentum fraction X 2 [^ . 
The PHENIX data are consistent with an increase of effective ‘absorption’ at forward rapidity, as discussed in Ref. Q, 
similar to that seen in fixed-target experiments at large Xp [5ll| . We will compare our results for shadowing alone 
with the PHENIX data at the end of this paper. 

Effects we do not consider here which affect the forward rapidity region in particular are energy loss in cold matter 
and intrinsic charm. The effect of energy loss, both with and without including shadowing, has been discussed in 
detail in Ref. (s^. Is^]. However, given the rather simple power law dependence employed for the p + p distribution, it 
would be worth pursuing embedding this approach within a production model such as the GEM. Intrinsic charm is 
not considered here because the interesting region, xp > 0.25, is inaccessible for quarkonium production in \y\ < 5, 
see Ref. 0. 

We do not include absorption by nucleons or comovers, which could affect the entire rapidity region, in this 
paper. The effective J/’ip absorption cross section has been seen to decrease with energy at midrapidity [^. When 
shadowing is included, the ^/s dependence is somewhat stronger due to the increased shadowing effect at low x [^ . 
If the nuclear crossing time is shorter than the J/4> formation time, the effective absorption decreases with 
as the state remains small during the entire time it spends in the target. However, absorption may play a role at 
the most backward rapidities when the quarkonium state is slow and can convert from the pre-resonance state to the 
final-state J/ip in the target [ 1 , [55|. The T(IS) absorption cross section is expected to be smaller than that of the 
J/V’ because of its smaller radius. 

While the J/'0-comover cross section is smaller than the nucleon absorption cross section, comovers may be more 
important for the excited charmonium and bottomonium states. (While these states have larger radii than the ground 
state quarkonium, this does not affect their potential absorption by nucleons since they still pass through the nuclear 
matter in their pre-resonant state.) Strong differences between the J/ip and ip' i?dAu at RHIC as a function of 
collision centrality would support the importance of comovers since their density increases with centrality. Recent 
comover-based p-|-Pb calculations show that this interpretation is consistent with the data [13] ■ 


V. NUCLEAR MODIFICATIONS OF THE PARTON DENSITIES 

We use several parameterizations of the nuclear modifications in the parton densities (nPDFs) to probe the possible 
range of shadowing effects: EPS09 [s^, EKS98 [il,l3j and nDSg [i^, and the two FGS sets, FGS-H and FGS-L 
[40l| . Since gg processes dominate quarkonium production over all measurable rapidities at the LHG, we focus more 
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on the behavior of the nuclear gluon parton densities in this section. Our main results in this paper are obtained with 
EPS09 NLO to make a fully consistent NLO calculation with the NLO CEM cross sections. However, in the following 
section, we will discuss the differences between the nPDF sets and between their LO and NLO manifestations on J/'0 
and T production in p+Pb collisions at = 5 TeV. 

All these sets involve some data fitting, typically nuclear deep-inelastic scattering (nDIS) data with additional 
constraints from other observables such as Drell-Yan dimuon production. The fact that these sets also include Drell- 
Yan production and do not include initial-state energy loss in matter in their global analyses, they exclude the 
possibility of additional quark energy loss since incorporating both over counts the effect on the sea quark densities. 
(Of course, if energy loss is generated by gluon emission either from an initial-state gluon or from the produced QQ 
pair, there is no energy loss possible in the Drell-Yan process, qq —>■ , since a virtual photon generates a lepton 

pair.) None of these data provide any direct constraint on the nuclear gluon density. It is thus obtained through fits 
to the dependence of the nuclear structure function, F^, and the momentum sum rule. The useful perturbative 
range of the nDIS data is rather limited since these data are only available at fixed-target energies. Thus the reach in 
momentum fraction, x, is also limited and there is little available data for x < 10~^ at perturbative values of p^. This 
situation is likely not to improve until an eA collider is constructed . We discuss the various nuclear modifications 
in order of their release. 

The EKS98 parameterization, by Eskola and collaborators [i^, IdJ, available for A > 2, is a leading order fit 
using the GRV LO 4^ proton parton densities as a baseline. The first EKS98 fits [i^ used the GRV LO proton 
densities from 1992 4^. This set employed a starting scale of Pq = 0.3 GeV. The set eventually released as EKS98 
was constructed with CTEQ4L [4l[ with initial scale of po = 1.6 GeV. To match the starting scale of the EKS98 
nPDF set, the CTEQ4L distributions were backward evolved from p = 1.6 GeV to 1.5 GeV They checked that 
the nPDF results from these two proton PDFs were consistent even though the proton PDFs themselves were quite 
different: the GRV LO set employed a valence-like gluon distribution at the starting scale, xg{x,fj.Q) ^ x°‘P{x) where 
a > 0 and P{x) is a polynomial function, while CTEQ4L takes a < 0, giving xg{x,fiQ) a finite value as a: —>■ 0. The 
minimum scales of the two sets are also very different. The kinematic range of EKS98, the only leading order set we 
consider since no NLO set was obtained at the time, is 1.5 < p < 100 GeV and 10“® < a; < 1. All the sets produced 
by Eskola and collaborators are based on piecewise functions at the minimum scale of the set that include small x 
shadowing, antishadowing at intermediate x, an EMC region (named for the effect first identified by the European 
Muon Collaboration) with < 1 at larger x and Fermi motion as x — >■ 1. 

deFlorian and Sassot produced the nDS and nDSg parameterizations at both leading and next-to-leading order 
for 4 < A < 208. Their results were based on the GRV98 proton PDFs at LO and NLO [^, using the same starting 
scales as GRV98, Pqlo = GeV, and Pqnlo = 0.25 GeV. GRV98 assumed valence-like input for the sea quarks 
and gluons at the minimum scales with xc(x,p§) = 0. The GRV98 gluon densities are not significantly different at 
low X when evolved to higher scales. The leading order set has a somewhat higher gluon density at low x [i^. Of 
the two deFlorian and Sassot sets, nDS and nDSg, the first, nDS, is an unconstrained fit that gives the best fit, 

To provide a set with stronger gluon shadowing, they constrained a second fit to give 5'^“ = 0.75 at x = 0.001 and 
p^ = 5 GeV^. This set, nDSg, was a much poorer fit to the New Muon Collaboration (NMC) and SLAG nDIS data 
and the E772 Drell-Yan cross sections [4^ . 

They showed that their results for the suppression factor for 7r° production at RHIC, 7?dAu(7'T), at both LO and 
NLO were in agreement with each other . This should be the case for a consistent order-by-order extraction of the 
nPDFs: not that the shadowing ratios are similar but that they give similar observable results, such as the nuclear 
suppression factors. The kinematic reach in x is the same as EKS98 while the p^ range is larger, 1 < p^ < 10® GeV^. 

The FGS sets by Frankfurt, Guzey and Strikman take the highest initial scale of all the nPDF sets, p = 2 GeV 
[40l| . They only provide NLO sets. Instead of a general global analysis, their nPDFs are based on the leading-twist 
approximation and rely on the diffractive PDFs. The Gribov-Glauber approach employed naturally introduces a 
centrality dependence for the shadowing. Their baseline proton parton densities are the CTEQ5M distributions. 
Therefore, their extraction is only good for x > 10“®. 

They extracted two sets to represent an upper and a lower bound on the shadowing ratios. FGS-H is gives stronger 
shadowing based on the two-nucleon scattering cross section. FGS-L is based instead on the ttN scattering cross 
section. Since this is larger than the two-nucleon scattering cross section, the shadowing effect derived in this case is 
weaker. Their approach is only applicable for x < 0.1. For x > 0.1, the sea quark ratio relative to the nucleon is fixed 
to unity and gluon antishadowing is obtained through the momentum sum rule. Hence, in this high x region, the 
two parameterizations are identical. Their approach does not allow predictions for shadowing on the valence quark 
distributions, these are taken from the LO EKS98 set. The FGS sets are available only for some specific nuclei where 
diffractive data exist: ^^G, ^®Ca, ^^®Pd, ^®^Au and ^®®Pb. The sets are valid in 10“® < x < 0.95 and 2 < p < 100 
GeV. 

The EPS09 [3^ LO and NLO parameterizations included uncertainties on the global analyses, both at LO and 
NLO, by varying one of the 15 fit parameters within its extremes while holding the others fixed. The upper and 
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FIG. 6: (Color online) Gluon shadowing ratios calculated for Pb nuclei {A = 208) calculated at the central value of the fitted 
factorization scales for J/psi (a) and T (b). The EPS09 NLO set is shown by the black solid curve while the uncertainty band 
is outlined by the black dotted curves. The NLO nDS and nDSg parameterizations are given in the blue dashed and dot-dashed 
curves respectively. The LO EKS98 parameterization is given in the magenta dot-dot-dash-dashed curve. The NLO FGS-H 
and FGS-L results are given by the red dot-dash-dash-dashed and dot-dot-dot-dashed curves respectively. 


lower bounds on EPS09 shadowing are obtained by adding the resulting uncertainties in quadrature [38|. The EPS09 


central LO results are in quite good agreement with the older EKS98 parameterization. The EPS09 LO sets are based 
on the CTEQ61L densities which are nearly flat at the initial scale as a:: —>■ 0. The EPS09 NLO densities are based on 
the CTEQ6M densities. Starting with CTEQ6M, a valence-like shape for the gluon distribution was adopted at the 
starting scale, in this case the charm mass threshold of p-o = 1.3 GeV [s^. The EPS09 sets assume the same starting 
scale as the CTEQ6 sets and is valid in the range 10“® < a; < 1 [s^. 

Figure El shows the sets we have discussed in this section at the scales of J/ip (a) and T (b) production. The results 
are shown for x values as low as 10“®. The EPS09 NLO sets are shown in the black solid and dotted curves. The solid 
curves show the central results while the dotted curves outline the NLO bands. The EPS09 results for J/xj} show a 
change in curvature, almost a dip, at x ^ 0.01, with a small rise at lower x followed by an eventual decrease. There is 
a rather strong antishadowing peak from 0.02 < x < 0.2 where the ratio drops below unity again, in the EMC region. 

The EKS98 LO ratio is similar to that of the central EPS09 LO set but decreases more smoothly, giving a stronger 
shadowing at lower x than the NLO set, a 40% shadowing effect at a; = 10“® rather than the 30% effect at NLO. At 
the lowest x values, the central EPS09 LO is equivalent to the lower limit of the EPS09 NLO band. 

The FGS-H and FGS-L sets are similar to the others but are somewhat narrower in the antishadowing region. 
However, they decrease faster with decreasing x than any of the other sets shown. They drop sharply at x = 10“®. 
Instead of either fixing the ratio at its value at a; = 10“®, as in the CTEQ5M proton PDF set, or trying to make a 
smooth extrapolation to lower x, it simply falls off more like a step function. 

The nDS set has somewhat stronger shadowing that the upper limit of the EPS09 NLO band while the nDSg set is 
between the EPS09 NLO central and lower limits. These sets, of all the results shown, have no antishadowing at all. 

The same features can be observed for the shadowing ratios at the T scale. The EPS09 band now decreases more 
smoothly as x decreases. The FGS-H ad FGS-L ratios show the steep drop at the limit where x = 10“® more clearly 
at the larger scale. All the slopes of the other nPDF sets become rather similar. 


VI. RESULTS 

We now show calculations of the J/'i/' and T production ratios, RpPh and Rfb as a function of rapidity and 
Pt, taking only nuclear shadowing effects into account. We begin with comparing results with the same shadowing 
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FIG. 7: (Color online) The J/'ip ratio RpPbiPr) for ALICE at forward rapidity 01 (a) and pr-integrated as a function of rapidity 
0 (b). The ratios are shown for CTIO (black), CTEQ6M (blue), CTEQ5M (red) and MSTW (magenta). 


parameterization calculated with different proton PDFs in Sec. IVI Al These results are shown for the J/^/> only and 
compared to a subset of the available data. In the remainder of this section, we compare the calculations to the full 
J/ijj and T data sets from ALICE and LHCb. We show the measurements compared to the EPS09 NLO uncertainty 
bands in Sec. IVI Bl Next, we compare and contrast the leading and next-to-leading order shadowing results of EPS09 
and nDS in Sec. IVI Cl Thereafter, the results from the nPDF sets shown in Fig. [S] are compared to the data and each 
other in Sec. IVI PI The mass and scale uncertainties on the central EPS09 set are compared to the uncertainties in 
the shadowing parameterizations themselves in Sec. IVI El Next we discuss how well the A + A results factorize into a 
product of p -I- A and A-\- p collisions in Sec. IVI FI Finally, we present the EPS09 NLO uncertainty bands for RHIC 
kinematics in Sec. IVI Gl 


A. Comparison of EPS09 NLO for Different Proton PDFs 

To begin, we check whether or not the shadowing results depend on the chosen set of proton parton distributions. 
As demonstrated in Sec. IIII Bl the choice of proton PDF has a strong effect on the shape of the individual p + p 
distributions. We now want to find out of this also means a change in the shape of the nuclear suppression factors 
since each shadowing parameterization is based on different proton PDFs. We employ the central EPS09 NLO set for 
this study. 

The results for RpPhipr) and RpPhiu) are shown in Fig. [T] The left hand side shows RpPhiPr) for the ALICE 
forward rapidity bin, 2.03 < j/cms < 3.53 while the right-hand side displays the pr-integrated RpPhiv)- The only 
visible difference between the results is for the lowest pT bin shown, 1 < p^ < 2 GeV, where the difference is ^ 4% 
between the ratios with CTEQ5M and MSTW relative to CTIO. There is no difference in the CTIO and CTEQ6M 
ratios. At higher px, the difference is negligible. Some variation on the order of the percent level is seen as a function 
of rapidity, particularly in the backward region. Since these variations are significantly less than those between the 
nPDFs themselves, we see that the choice of proton PDF has a negligible effect on the results. 


B. EPS09 NLO Uncertainties 

Here we present the EPS09 NLO uncertainties in J/ ip and T production. We will compare the calculated nuclear 
suppression factor, i?ppb, and the forward-backward ratio, Rfb, to the ALICE and LHCb data. 
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FIG. 8: (Color online) The J/%p ratio RpPbipr) in the ALICE acceptance at forward (a), backward (b) and central (c) rapidity 
Ql- The ratio RpPb{y) [H is shown in (d). The EPS09 NLO uncertainty band is shown. 


1. J/V> 

Figure H] compares the suppression factors for J/ip ^ as a function of pr at forward (a) and backward (b) 

rapidity as well as for J/?/; —> e“'"e“ at midrapidity (c). The rapidity dependent ratio is shown in (d). The solid curves 
show the EPS09 NLO central results while the dashed curves outline the upper and lower limits of the EPS09 NLO 
uncertainty band. In all cases, the bands on the EPS09 NLO uncertainties are obtained by adding the differences in 
the results obtained by varying each parameter separately by one standard deviation of the fit in quadrature. 

In general, the pr-dependent data are in relatively good agreement with the EPS09 NLO bands, if not with 
the central values themselves. In the forward region, the measured i?ppb is compatible with the lower edge of the 
uncertainty band for pT < 6 GeV. The data at backward rapidity also suggest a somewhat stronger cold matter effect, 
compatible with the upper edge of the band. While the uncertainties on the data at midrapidity are larger than 
those on the calculation, the data again indicate a somewhat stronger effect than suggested by EPS09 NLO. However, 
within the uncertainties of both the data and the calculations, the results are, in general agreement. 
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FIG. 9: (Color online) The J/ip forward-backward ratio Rfb{pt) in the ALICE overlap region, 2.96 < j/cms < 3.53 (a) and 
RpBiy) (b). The EPS09 NLO uncertainty band is shown with the ALICE [H and LHCb [j data. 


A similar conclusion can be drawn from the pr-integrated rapidity dependence. The ALICE data at forward and 
backward rapidity are shown in two ways: a single pr-integrated point over the entire rapidity acceptance and with 
each broad y bin broken up into six separate bins. While the smaller bins are almost independent of y in the backward 
region, there is a decrease in Rpphiu) in the forward region. The overall observed dependence of RpPhiv) is in good 
agreement with the pT-dependent results in Fig. [5Ka)-(c). While the LHCb are within one standard deviation of the 
ALICE data, they are somewhat lower. 

As discussed previously, the suppression factor Rpph is somewhat artificial since there is currently no p+p reference 
measurement at the same energy as the p-|-Pb data. Instead the reference is obtained from an interpolation between 
the p + p measurements at 2.76 and 7 TeV, as described in Sec. m The forward-backward ratio reduces the systematic 
uncertainty. 

The results are compared to the J/ijj calculations in Fig.[9l We note that the uncertainties on the ratio are narrower 
than those of the forward and backward regions separately, as shown in Fig. [51 The uncertainty bands are formed by 
taking the forward-backward ratios for each of the 31 EPS09 sets individually and then adding the uncertainties in 
quadrature. The calculated ratios are almost independent of pt- This behavior is due to the fact that the forward 
ratio is less than unity and increasing with px while the backward ratio is greater than unity over most of the pt 
range. The data instead show a minimum at low px, increasing to RpPh ^ 0.8 at pt ^ 8 GeV, albeit with large 
statistical uncertainties. 

The calculated ratios as a function of rapidity are very narrow for y < 2.5 before broadening at larger y. This 
behavior is obvious from looking at RpPh{y) since, for y > —2.5 the uncertainty band is parallel to the central value and 
almost linear so that when the forward-backward ratio is formed, the band is compressed. The ratio Rpsiy) broadens 
and decreases at p > 3, where, at backward p, RpPhiy) enters the antishadowing region. The forward-backward ratio 
narrows again at p ~ 4, near the ‘pinch’ in i?pPb(p), where the data are best constrained, and finally broadens again 
in the antishadowing regime where the constraints are poorer. The data are below the calculated band because, while 
the backward region is rather well described, the central value is above the data at forward rapidity. We note that 
while the ALICE RpB{y), in a smaller rapidity range, is almost flat, the LHCb ratio varies more. However, in the 
bin where the rapidity ranges of the two experiments overlap, they are in good agreement. 

Finally we remark on the recent ATLAS forward-backward measurement of Rpb{pt) in the central region, |pcms| < 
1.94 and RpBiy) for 8 < pt < 30 GeV [l^. Their results are consistent with unity in both pp and p and agree 
well with our EPS09 NLO calculations. The results shown in Fig. Ha) are in the forward region while the ATLAS 
measurement is at midrapidity. Given the higher scale implicit in the pp range covered by ATLAS, as well as the 
higher x in the rapidity range covered by ATLAS, it is unsurprising that the EPS09 NLO result is consistent with 
unity here. Similarly, given that the EPS09 NLO result integrated over pp in Fig.j^b) is already similar to unity for 
IpI < 1.5, we can expect that the same ratio for pp > 8 GeV has a more shallow slope and deviates less from unity 
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than the results compared to ALICE and LHCb here. 


2 . T 

The results for the T suppression factor, RpPh, are shown in Fig. (TUI In this case, the data are insufficient to form 
RpPhipr) or Rfb(jPt)- In addition, no midrapidity e+e” T measurement has been reported. It is also not possible 
to separate RpPhijj) into smaller bins due to the low statistics. Finally, the ALICE and LHCb T data, in the forward 
and backward rapidity bins, while within one standard deviation of each other, do not appear to be in very good 
agreement with each other. Therefore, it is more difficult to draw conclusions about the T results. 

The pT-dependent ratios differ from those for J/ij} most at low pT- In particular, the backward rapidity region shows 
an antishadowing effect already at low pT- At higher pT values, the scales are similar since pt > to in both cases. 
The EPS09 NLO uncertainty band is somewhat narrower, especially as a function of rapidity, where the ALICE and 
LHCb T data are both shown. The agreement with the LHCb data is quite good while the agreement with ALICE is 
rather poor. The two data sets are approximately within one standard deviation of each other but LHCb is indicative 
of antishadowing at the approximate rapidity of the EPS09 NLO antishadowing peak while the ALICE result is not. 
At forward rapidity, the ALICE J/ij) and T results for RpPh are similar. 

The T forward-backward ratio is shown in Fig. [11] The pr-dependent ratio is lower than that of J/tp due to the 
greater T antishadowing. The uncertainty is larger because the uncertainty on the backward rapidity J/ip RpPh is 
narrower than for T since the ALICE J/ip region is near the ‘pinch’ between the shadowing and antishadowing regions 
while the T rapidity is directly in the region where the uncertainties are larger for gluons. This is also reflected by 
the rapidity dependence. 


C. Comparison of Leading and Next-to-Leading Order Results 

Previously we reported that the shadowing parameterizations gave the same results for the LO and NLO cross 
sections [s^. That result was based, however, on employing the LO EKS98 set in the LO and NLO CEM calculations. 
In addition, other authors have suggested that the order of the calculation does not matter and NLO nPDF sets can 
be used with LO quarkonium calculations. It is worth checking these assumptions in detail. 

As discussed in the previous section, the nDS and nDSg LO and NLO sets were checked for consistency for the 
PHENIX 7r° data [i^. We now do the same for J/ip and T production as a function of rapidity for EPS09 LO and 
NLO and for nDS and nDSg LO and NLO in Fig. [T^ 

The shadowing parameterizations are compared in Fig.[T^(a) and (b). There are considerable differences between 
the EPS09 bands at a; < 0.01. The NLO result has a different curvature than the LO one, changing slope and 
becoming flatter, decreasing slowly as x is lowered. On the other hand the LO bands decrease smoothly. There is, 
however, very little difference in the upper limit of the band, with the weakest shadowing effect. The difference is 
more pronounced for the central sets and striking for the lower limit of the band, with the strongest shadowing. In 
this case, the lower NLO limit is close to the central EPS09 LO set. There is also some difference between the nDS and 
nDSg parameterizations. While the difference is essentially negligible for the nDS set, there is a stronger deviation 
between the nDSg curves, particularly around x ^ 0.01. For x < 0.002, the ratios are parallel. 

We now turn to whether these differences affect RpPh{y) for quarkonium production. The EPS09 LO calculation is 
made with a fully LO CEM calculation in 2 —>■ 1 kinematics. The NLO calculation is with the NLO matrix elements 
in the CEM. 

Results are shown for in Fig. [T^ (c) and (d). There is clearly a strong effect for EPS09 with the central LO 
result passing through the ALICE and LHCb data while only the lower limit of the NLO band reaches the data, an 
ap-fproximate 15% difference between the two results. However, the nDS and nDSg results, while clearly not on top 
of each other, as they are for the 7r° calculation [d^, agree within ^ 3%. 

Why is one derivation nearly consistent order by order and the other is not? We suspect that the primary cause is 
the low X behavior of the baseline proton parton densities, as already mentioned in Sec. [Vj The CTEQ6M, GRV98 
LO and GRV98 NLO gluon distributions are all based on valence-like initial distributions at po, x°‘P{x). However, 
the CTEQ61L gluon density takes an almost constant finite value at i.e. a ^ 0. This initial difference carries 
forward over all p^ because the NLO set must have stronger p^ evolution to fit the same sets of DIS data as at LO. 
Thus the CTEQ6M distribution evolves much more rapidly in p'^ than does CTEQ61L. At the factorization scale of T 
production, the difference between CTEQ6M and CTEQ61L is still substantial while the GRV98 LO and NLO gluon 
distributions, with similar starting behavior, are closer together. The EPS09 LO shadowing thus has to be stronger at 
low X to produce the same behavior at higher scales as the EPS09 NLO shadowing. The compensation does not have 
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Pt (GeV) 


FIG. 10: (Color online) The T ratio RpPh{pT) in the ALICE acceptance at forward (a), backward (b) and central (c) rapidity. 
The ratio RpPh{y) @| is shown in (d). The EPS09 NLO uncertainty band is shown. 


to be as large for the nDS sets since the GRV98 LO and NLO sets are more similar at low x and the scale relevant 
for quarkonium production. 

The agreement between the EPS09 LO and NLO results is not improved for the higher scale of T production, as 
shown in Fig. [12] (e) and (f). The higher scale, as already noted, weakens the shadowing and narrows the uncer¬ 
tainty bands but does not bring the central results at LO and NLO into better agreement. This conclusion is likely 
independent of the specific final state as long as its production is dominated by low x gluons. 

There are, however, valid reasons for there to be a few percent difference between the LO and NLO results seen 
with nDS and nDSg. The LO OEM calculation is a 2 —>■ 1 process with fixed values of Xi and x^ for a given y since 
there is no pt scale in the calculation, only in the factorization scale which can be adjusted to include an average 
value of Pt- The NLO calculation, on the other hand, includes 2 —>■ 2 and 2—^3 processes (the LO-I-virtual and real 
contributions respectively) to the NLO GEM result. The NLO GEM calculation does not have a fixed correspondence 
between xi, X 2 and y because of the pr scale. This higher scale can also lead, on average, to a larger factorization 
scale and a somewhat larger X 2 on average in the NLO calculation. All these differences can cumulatively lead to the 
2 — 3% effect between LO and NLO observed for the nDS sets. Similar arguments can also explain the differences 
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FIG. 11: (Color online) The T forward-backward ratio Rfb{pt) in the ALICE overlap region (a) and Rpsiy) (b). The EPS09 
NLO uncertainty band is shown, along with the ALICE and LHCb data. 


between the LO CEM and LO CSM results and are not directly attributable to the production mechanism per 
se. 


D. Comparison of Shadowing Parameterizations 

In this section, we compare the central EPS09 NLO results to those with other shadowing parameterizations. 


1. J/^ 

Figures [13] and |T4| compare the NLO nDS, nDSg, FGS-H and FGS-L results to the central EPS09 NLO result 
(shown in black). The results with the EKS98 LO set are also shown. In the pr ratios, the EPS09 central set tends 
to underestimate the effect relative to the data except at backward rapidity. Indeed, at forward and midrapidity, the 
nDS set is the only one that underestimates the data more since it has the weakest gluon shadowing. At backward 
rapidity, the nDS result is close to unity over all pt while nDSg exhibits some shadowing. This is because these two 
sets have no gluon antishadowing. At forward rapidity, the FGS-H and FGS-L sets result in the strongest effect, 
significantly stronger than the data in the case of FGS-H. The set that comes closest to agreement with the data in 
all three rapidity regions is EKS98 which only overestimates the shadowing effect for px > Q GeV. 

In all cases the pT-dependent results for all the nPDF sets shown agree within 10% for px > 10 GeV. This can be 
expected because the evolution of the modifications is relatively strong for the gluons. 

As seen in Fig. lI3f a). the spread between predictions is largest in the forward region where the shadowing predictions 
differ most. Indeed, in the px bin centered at 1.5 GeV, there is a factor of ^ 8 between the values of RpPh- The gap is 
reduced to ^ 1.3 by px ^ 3.5 GeV. The weakest px dependence is given by the calculations with nDS and the central 
EPS09 NLO set. The nDSg set also results in a weak px dependence but has a stronger overall shadowing effect. The 
strongest shadowing comes from the FGS sets which overpredict the shadowing strength. The best agreement with 
the ALICE data is obtained for the EKS98 LO set which has no NLO counterpart. We note that the measured RpPh 
approaches unity a,t px 7 GeV while none of the shadowing parameterizations give a result approaching unity over 
the entire px range shown. 

In the backward region, illustrated in Fig. [T5I b). the results show some antishadowing, except for nDS and nDSg 
which specifically exclude it. All the calculations lie below the centroids of the data, partly because the ALICE 
acceptance is centered on the side of the antishadowing peak so that the calculated antishadowing is, on average. 
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FIG. 12: (Color online) The EPS09 central results, as well as the uncertainty bands, are shown in (a), (c) and (e). In all 
cases, the solid red curve shows the central NLO result while the dashed red curves delineate the NLO uncertainty band. The 
dot-dashed blue curve is the LO central result while the dotted curves outline the LO uncertainty band. The nDS and nDSg 
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Pt (GeV) 


FIG. 13: (Color online) The J/'ip ratio RpPhipr) in the ALICE acceptance at forward (a), backward (b) and central (c) 
rapidity Q|. The ratio RpPh{y) for ALICE [II] and LHCb 0] is shown in (d). The results are shown for central EPS09 NLO 
(black), nDS NLO (blue dashed), nDSg NLO (blue dot-dashed) and EKS98 LO (magenta dot-dot-dash-dashed), EGS-H (red 
dot-dash-dash-dashed) and FGS-L (red dot-dot-dot-dash). 


lower than that of the data. Note also that the antishadowing peak is broad in px, with some antishadowing remaining 
at Pt 15 GeV. 

The pT-dependent results from the nPDF sets come closest together in the midrapidity bin, with a spread of at most 
30%, as seen in Fig. IT^ c). Although the data again suggest stronger shadowing at midrapidity than the calculations, 
the uncertainties in the data are large since the number of J/ip in this bin are lower by more than a factor of 100 
than in the forward and backward bins, see Table ID In this rapidity range, the shadowing effect is reduced to less 
than 10% for px > 10 GeV with all sets. 

The largest difference in the RpPh ratios is as a function of rapidity, shown in Fig. [T3l d). since the pr-integrated 
ratios access the lowest x values. The ratios are closest at backward rapidity, not surprisingly, because the EPS09 
NLO and EKS98 sets follow each other closely for x > 0.002, see Fig.jSKa). Since the FGS sets are somewhat narrower 
in the antishadowing x region, this is manifested as an apparent shift toward negative rapidity in the antishadowing 
peak for the J/ip. The EPS09 NLO and EKS98 are in closest agreement with the integrated ALIGE rapidity bin which 
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FIG. 14: The J/'i/) forward-backward ratio Rfb{pt) in the ALICE overlap region (a) and Rpsiy) (b) [J]. The ratios are for 
central EPS09 NLO (black), nDS NLO (blue dashed), nDSg NLO (blue dot-dashed) and EKS98 LO (magenta dot-dot-dash- 
dashed), FGS-H (red dot-dash-dash-dashed) and FGS-L (red dot-dot-dot-dash). 


is larger than unity. The ALICE data were split into six smaller bins in both the forward and backward rapidity bins. 
In the backward region, the smaller bins are almost independent of rapidity which is not consistent with any of the 
shadowing results. Since the LHCb backward rapidity result is below unity, only the nDS result is not in agreement 
with LHCb. 

The ratios separate further in the shadowing regions at mid and forward rapidity. Since the ALICE pr-integrated 
ratio at midrapidity is almost the same as that at forward rapidity, albeit with somewhat larger uncertainty, the 
midrapidity point is only in agreement with the FGS-L and FGS-H sets. This is not surprising: because of their steep 
drop at low x, they show stronger shadowing than EKS98 already at a: ^ 0.001 and considerably overestimate the 
effect at forward rapidity where x < 10“^. In the forward bin, the EKS98 and nDSg results are in best agreement 
with the ALICE and LHCb results. The ALICE forward bin, when split into smaller bins, shows a linear dependence 
on rapidity but with a slope steeper than any of the shadowing calculations. Overall the EKS98 parameterization, one 
of the oldest and LO only, agrees best with the data. Since it is similar in shape to the EPS09 LO central set, shown 
in Fig. IT^ a). this is not surprising. The EPS09 LO and EKS98 sets, based on CTEQ6IL and CTEQ4L respectively, 
are the only sets of those from EPS09 and nDS that are not based on valence-like gluon distributions in the proton 
and are thus subject to slower scale evolution. 

In Fig. [Ml the forward-backward ratios are shown as functions of pr and y. The EPS09, nDS and nDSg ratios are 
nearly independent of pr, thus above the data over most of the pr range. The EKS98 and FGS sets have a somewhat 
stronger dependence on px and thus agree more closely with the forward-backward ratio at low px- However, of 
these, only the FGS-L and EKS98 sets are in agreement with the rapidity-dependent ratio. All the sets, with the 
exception of the central EPS09 NLO set, give forward-backward ratios that are linear in y. As already discussed this 
is due to the abrupt change in slope of EPS09 NLO at x ^ 0.002, see Fig. UDa), at the transition from shadowing to 
antishadowing, absent from the other sets shown. 

Overall, the EKS98 LO set seems to agree best with the J/ij} data, both i?ppb and Rfb- This is somewhat less 
than satisfying since we are using a LO nPDF set with a NLO cross section calculation. In addition, it is one of the 
older sets employed here. However, it agrees rather well with EPS09 LO when calculated with LO matrix elements, 
even though the baseline proton PDF and its behavior at low x is very different. We had previously demonstrated 
that the ratio RpA{y), calculated at LO and NLO with the EKS98 LO set, gave almost identical results at each order 
[s^ . This is not a surprise since the gg contribution dominates production both at LO and NLO and the shift in x 
and due to the different scales is small when integrated over px- 
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FIG. 15: (Color online) The T ratio RpPh{pT) in the ALICE acceptance at forward (a), backward (b) and central (c) rapidity. 
The ratio RpPh{y) is compared to the ALICE and LHCb data in (d). The ratios are for central EPS09 NLO (black), 
nDS NLO (blue dashed), nDSg NLO (blue dot-dashed) and EKS98 LO (magenta dot-dot-dash-dashed), FGS-H (red dot-dash- 
dash-dashed) and FGS-L (red dot-dot-dot-dash). 


2. T 

Figures Hn] and Uni present the results for T production in the same rapidity ranges. The pT-dependent Rppp, ratios 
shown in Fig. m a)-(c) exhibit somewhat less spread than the counterpart J/ij} calculations. The main difference is 
that now, at the higher scale and correspondingly larger x probed, the antishadowing peak for T production is fully 
within the acceptance of ALICE and LHCb in the backward direction. Thus RpPhiPr) > 1 over all px in Fig. [TSl b) 
with RpPhipr ^ 1.5 GeV) ~ 1.2. 

The rapidity-dependent ratio, RpPh{y), shown in Fig. [T5T d) echoes the low px results. Even the nDS and nDSg 
ratios are above unity in the backward rapidity region, albeit in the backward edge of the experimental acceptance. 
This is because, in the high-a: region, these sets become large as a: —^ 1 with the turn on of this effect beginning at 
X ^ 0.5 and the region of x ^ 0.1 is where the gluon distribution becomes larger than unity. We note that even 
though the FGS gluon ratios are only somewhat narrower than the EPS09 NLO ratio for the same scale, see Fig.jSKb), 
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FIG. 16: The T forward-backward ratio Rfb{pt) in the ALICE overlap region (a) and Rfb{v) (b). The ALICE and 
LHCb data are shown in (b). The ratios are shown for central EPS09 NLO (black), nDS NLO (blue dashed), nDSg 
NLO (blue dot-dashed) and EKS98 LO (magenta dot-dot-dash-dashed), FGS-H (red dot-dash-dash-dashed) and FGS-L (red 
dot-dot-dot-dash). 


for the actual T calculation shown here, the difference is enhanced. At backward rapidity, the sets are in agreement 
only with the higher LHCb point while they all miss the ALICE point at RpPh{y) <1. At forward rapidity, the 
results spread more. The strongest shadowing is again with the FGS results but due to the larger x probed for the 
T, the X < 10“^ region is not entered in the rapidity range shown although it is approaching the edge of this region, 
contributing to the fluctuations in the FGS curves in Fig. ITSf d). As for the J/ip, the EKS98 slope is more linear for 
y > —1 than than the EPS09 NLO central value. Thus the nDS and nDSg results, with weak shadowing, as well 
as the EPS09 NLO results agree best with the weaker shadowing reported by LHGb while the FGS results generally 
agree better with the stronger shadowing reported by ALIGE. The EKS98 LO parameterization is between the two 
measurements and also includes antishadowing. Therefore, it gives the best overall result for this measurement given 
the quality of the currently available data. 

The forward-backward ratios for T production are shown in Fig. 1161 None of the results for the various nPDF sets 
exhibit a strong pT dependence. The FGS and EKS98 results are generally below those of EPS09 NLO over all px 
while those of nDS and nDSg are higher. These rather flat ratios can be attributed to the relative pT-independence 
for pt > i GeV for RpPhipr) in Fig. fW a) and (b) while the stronger shadowing and antishadowing, respectively, in 
these regions at pt < 4 GeV essentially give the same forward-backward ratio as at higher px- 

The rapidity dependence of the forward-backward ratio. Fig. ITST b). is similar for all the nPDF sets even though 
the magnitude of the effect varies. Because all the ratios increase above unity in the backward direction, the ratios 
all decrease almost linearly, except for EPS09 NLO, with a minimum at 3.5 < y < 4.0. While most of the ratios rise 
again at larger rapidity, the nDS and nDSg ratios only flatten with rapidity above y ^ 3.5. 

In summary, of all the nPDFs, the LO EKS98 set agrees best with the currently available data. However, this is a 
rather unsatisfying conclusion because the LO set is used in a NLO calculation. Even if the x values do not strongly 
differ between LO and NLO, the LO shadowing parameterization should be a stronger function of x, with a larger 
overall shadowing magnitude, to produce the same RpPh ratios order by order, as was previously shown to be the case 
for the nDS and nDSg sets. A re-evaluation of the nPDF analysis, including appropriate LHC data, should be made 
in order to resolve the issue. 


E. Mass and Scale Dependence 

In this section, we discuss the mass and scale dependence of the shadowing ratios, using the EPS09 NLO 
set as an example. Similar results are expected for the other nPDFs but we choose EPS09 NLO here be- 
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cause we can compare the shadowing uncertainty to that of the mass and scale. Recall that for charm- 
we employ {m, fip/m, fj-fi/m) = (1.27 ± 0.09 GeV, 2.lto;85i l-6lo;i2) while for bottom-T production, we employ 
(m, fip/m, irnlm) = (4.65 ± 0.09 GeV, IA^^oaI l-ltoJo)- 

The uncertainties of the EPS09 NLO nPDF sets compared to those of the mass and scale uncertainties on da/dpT 
for the J/ip are illustrated in Fig. [T71 Here the ALICE pt distributions in p-|-Pb collisions Q are compared to our 
NLO GEM calculations. The red solid and dashed curves are the central EPS09 NLO results with the uncertainty 
band due to the EPS09 parameter variations. This band is rather narrow on the scale of the plots with the only 
clear separation between the curves at low pT where the shadowing effects are largest. Clearly, a 20-40% effect on 
RpPbipr) on a linear scale in this region results in a narrow band on the px distribution itself on a logarithmic scale. 
We do not show the lowest pr results because the un derly ing p + p calculation does not accurately numerically cancel 
the divergences in the negative weight Monte Carlo [^. The intrinsic kr kick used to obtain the shape of the pr 
distribution in the low-px region was fixed at RHIC energies and successfully compared to the J/ij} px distributions 
at = 7 TeV [H]. 

The magenta curves in Fig. [T7] show the uncertainty due to the mass and scale variations. The mass and scale 
uncertainties are calculated based on results using the one standard deviation uncertainties on the quark mass and 
scale parameters. If the central, upper and lower limits of are denoted as C, H, and L respectively, then 

the seven sets corresponding to the scale uncertainty are {{fip/'m, dp/'m)} = {{C,C), {H,H), {L,L), {C,L), {L,C), 
{C, H), {H, C)}. The uncertainty band can be obtained for the best fit sets by adding the uncertainties from the mass 
and scale variations in quadrature. The envelope containing the resulting curves, 

dUniax/dA = dfJcent/dA -\~ ^ (dfJ^^niax/dA dfJcent/ {d(Jm ,maK /dJC dfJcent/dA)^ , (5) 

du'min/ dX = dfJcent/dA ^ (dcT^^min/ dX do'cent/dX')‘^ -j- {d<7m,inin/ dAl dcT^eiA/dX^'^ , (6) 

defines the uncertainty on the cross section where X can be either px or y. The EPS09 uncertainty band is based on 
the mass and scale set with the central mass value and {pF/fn, dp/ni) = {C,C). 

In Fig. 1171 the mass and scale uncertainties are clearly significantly larger than those of the shadowing parame- 
terizations alone, still a factor of ^ 2 on the J/ij} px distributions over a wide px range. Thus, even though the 
uncertainties on the J/'0 cross section in the GEM are reduced relative to those calculated previously, see Ref. [l^ . 
they are still significantly larger than those of shadowing on the central mass and scale values. One would naively 
expect the uncertainties on i?ppb and Rfb to reflect the results shown in Fig. [171 but, as we will show, it depends on 
the method used to calculate the uncertainty. 

We have tried several ways of estimating the size of the mass and scale uncertainty on the ratios which we describe 
first before comparing to the data. We note that all of these calculations are based on the EPS09 NLO central set 
rather than making a global nPDF -|- mass -I- scale uncertainty by calculating the 30 error sets for EPS09 with the 8 
additional mass and scale combinations that define the mass and scale uncertainty on the cross section for p p and 
p-|-Pb. 

The most naive, labeled m/ pf/PrvI on Figs. ITBIITTI is to take the ratios of p-|-Pb top-|-pfor each mass and scale 
combination and then locate the extrema of the ratios based on the mass and scale in each case, exactly as in Eqs. 
and dni above. Thus the uncertainty band is based on the RpPh of each set. Since the ratios are all close to unity 
(within 20-40% in most cases), and similar to each other, this method would tend to underestimate the uncertainty. 

The second, labeled m/ pf/Prv2 on Figs. [TB1I771 reflects the uncertainty calculation of Eqs. ([S]) and dH). We 
calculate the maxima and minima of da/dX, as in Fig. 1171 and then form Rppb by dividing by the central value of the 
p+p cross section in the corresponding rapidity bin (for the pT-dependent results) or the pp-integrated cross section 
(for the rapidity dependence). Thus the shadowing ratios for v2 are based on cross sections only. We expect a larger 
uncertainty band for this calculation, especially at low px- 

We note, however, that there is actually no difference between the v2 and z;l calculations of the uncertainty in 
the forward-backward ratios because we calculate the forward-backward ratios for each mass and scale combination 
before calculating the uncertainty to ensure that we make a one-to-one correspondence in Rfb- Thus, there is no v2 
result in the plots of Rfb- 

The above methods of calculating the mass and scale uncertainty, relying on Eqs. and ®, are not actually 
calculated the same way as the EPS09 uncertainty band because they rely only on the extrema of the mass and 
scale dependent cross sections from the central values rather than adding the excursions from the central value in 
quadrature, as in the EPS09 uncertainty calculation. Therefore, finally, for m/ pf/Prv3, we add the mass and scale 
uncertainties in quadrature, a la EPS09, and then form RpPh by dividing by the central p -|-p cross section. Since this 
is a cumulative uncertainty rather than based on the greatest excursion from the mean, we expect this to give the 
largest overall uncertainty. The v3 uncertainty band was calculated assuming that the appropriate pf/tti and pR/m 
pairs are [{H, H), (L, L)], [{H, C), {L, C)] and [(C, iJ), (C, L)]. Other choices could lead to different results. (There is 
only one possible pair for the mass uncertainty thus this part of the calculation is the same for v2 and v3.) 
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FIG. 17: (Color online) The ALICE J/ip pr distributions [2] at forward (a), backward (b) and midrapidity (c) at NLO in the 
CEM [T^. The solid red curve is the EPS09 NLO central value while the dashed red curves are the EPS09 NLO uncertainties 
and the dot-dashed magenta curves are the mass and scale uncertainties. 


Results using all three methods are shown for J/ip and T in the remainder of this section. 


1. J/ip 

Figure [T5] shows the relative uncertainties for the EPS09 NLO band and the three ways of calculating the mass and 
scale uncertainty for RpPh- As expected, the t;l proceedure gives the smallest uncertainty since there is little variation 
in the individual values of RpPh for each mass and scale choice for the central EPS09 set. Therefore, the v\ band is 
narrower than that of the EPS09 band itself, underestimating the uncertainty. 

On the other hand, the v2 and v2> bands are wider, especially for pp < 5 GeV, as expected. Also, as mentioned 
previously, the vS band is broader than that of v2 over all pp although the two methods merge at the highest pt 
values. The RpPhiPx) found for the v2 method is in agreement with what one might expect looking at the bands 
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FIG. 18: (Color online) The mass and scale uncertainties in the J/'ip ratio RpPhijpr) in the ALICE acceptance at forward (a), 
backward (b) and central (c) rapidity Q|. The ratio RpPh{y) is compared to the ALICE [I[ and LHCb 3 data in (d). The 
EPS09 NLO uncertainty band is shown in red while the uncertainties calculated with method ul in blue, v2 in magenta and 
u3 in black. 


on da/dpT in Fig. [T71 The ratios as a function of pr are relatively regular even though, at least for v\ and u2, the 
combination of factorization and renormalization scales that give the extreme values for that pt does depend on pT- 
The situation is somewhat different as a function of rapidity, in part because the rapidity distributions fluctuate 
somewhat, especially at midrapidity. The size of the fluctuations depends on the scale choice and, depending on 
which scale set determines the extreme value, these fluctuations can manifest themselves in Rpp\j{y). This effect is 
obvious in Fig. IlSf d) where the upper bound on vl is visibly above the central value of the ratio for \y\ > 2 and 
almost on top of it for \y\ < 2. The abrupt change of slope occurs because the maximum value of the ratio is obtained 
with set {H,H) away from midrapidity and with set {C,L) at midrapidity. This is illustrated in Fig. [TH] where the 
p + p and p+Pb rapidity distributions are shown for each of these combinations. The change in the p+Pb rapidity 
distribution is a stronger function of y for the (iL, H) set than for the (C, L) set. There are also larger fluctuations in 
the distribution at midrapidity for the {H,H) set. To emphasize this difference, we present the rapidity distributions 
as histograms rather than smoothed curves. We note that these fluctuations are not due to limited statistics in the 
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FIG. 19: (Color online) The p + p and p+Pb J/ijj rapidity distributions for the {H,H) {C,L) sets showing the differences 
leading to the change in the upper limit of the mass and scale uncertainties of method vl around midrapidity. 


calculation but are more likely due to the uneven distribution of low pT negative-weight events in the MNR code at 
midrapidity. The presence of these events may also flatten the rapidity distribution in the LHC energy regime so that 
while the p + p rapidity distribution is within the mass and scale uncertainty band the slope is somewhat flatter than 
the measured one, see Ref. [l^ . 

The v2 and v3 ratios are more separated and, because they are obtained by taking the ratio to the p+p distribution 
calculated with the central mass and scale set, they are also smoother. The ‘kink’ in the central EPS09 NLO ratio 
noted previously is sharpened for the upper limit of these bands and reduced, but still present, in the lower limits. 
The upper limits obtained for v2 and u3 are very similar. Indeed these upper limits are almost indistinguishable for 
most of the illustrated rapidity range. The numerical values are generally different, with the u3 values larger than 
those of v2, but the difference is not large enough to be visible on the scale of the plot. This similarity likely stems 
from the aforementioned fluctuations in the rapidity distributions. Since these fluctuations tend to be smaller for the 
lower limits of Rpphiv), there is a large separation between the lower bounds on v2 and u3. 

In all the cases shown, except for the ratios calculated with vl, the envelope of ratios described by the mass and 
scale uncertainties contains the data. This can be expected because this is the designed purpose of Eqs. m and m 
for the p + p distributions and, when applied correctly, it does the same for the nuclear suppression factor Rpph as 
well. 

This is not necessarily the case for the forward-backward ratios, shown in Fig. Indeed, in this case the ratios 
calculated with methods ul and u3 are nearly identical. Recall that there is no difference between ul and v2 uncertain¬ 
ties for the forward-backward ratios. The most interesting thing to note here is that the ratios of the mass and scale 
uncertainty as a function of rapidity have different slopes than those of the EPS09 NLO uncertainty. The stronger 
rise in RpPhiy) at backward rapidity, combined with the weak dependence on rapidity in the forward direction makes 
the upper limit on the mass and scale dependence larger than the central EPS09 set at y < 2. At larger rapidities, 
these values are almost coincident. On the other hand, the weaker rise in RpPhiy) for the lower limit in the backward 
direction makes the lower limit of the forward-backward ratio decrease more rapidly than the central set for y > 3. 


2. T 

The results for the mass and scale uncertainties on the ratios for T production are presented in Figs. [21] and 
\n[ Overall the trends are quite similar to those for the J/'0 even though the pT-dependent results here are almost 
independent of pt- However, we note that the results for Rpph show much smaller differences between v2 and u3 over 
the entire px range, as well as a function of rapidity. In this case also, the shape of the limits on the mass and scale 
uncertainty bands for RpPh{y) is the same as that of the central EPS09 NLO value so the shape of that ratio has a 
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FIG. 20: (Color online) The mass and scale uncertainties in the J/'i/) forward-backward ratio Rfb{pt) in the ALICE overlap 
region (a) and Rpaiv) (b). The EPS09 NLO uncertainty band is shown in red while the uncertainties calculated with method 
vl in blue and u3 in black. The ALICE [J| and LHCb Q] data are also shown. 


weaker scale dependence than the J/ip- The band on RpPhiv) is broad enough to encompass the range of the low 
statistics T data. 

The uncertainties on the forward-backward ratio in Fig. [22] are again small. They are, in fact, considerably smaller 
than those due to EPS09 which, as we saw previously, is not because the uncertainties on the distributions in an 
individual rapidity bin, for Rfb{pt), are small. Instead, these uncertainties are somewhat canceled in the ratio, even 
though we still expect the largest excursion from the central value for u3. 


F. Factorization 

The question of whether cold nuclear matter effects factorize has been much discussed. If so, then at the same 
energy, the product of the nuclear modification factors at forward and backward rapidity in pA collisions would give 
the predicted cold matter result for A + A collisions at the same energy. 

In their latest paper on the J/ij} p-|-Pb results 0, the ALICE Collaboration compared the product of their RpPhiPr) 
ratios in the forward and backward rapidity regions at ^/s^ = 5.02 TeV, 2.03 < j/cms < 3.53 and —4.46 < Ucms < 
—2.96 respectively, to the i?pbPb(PT) ratio in the region 2.5 < y < 4.0 (symmetric around midrapidity) at = 2.76 

TeV. They did the same for the square of the midrapidity ratio measured in —1.37 < ?/cms < —0.43 in p-|-Pb collisions 
= 5.02 TeV to i?pbPb in \y\ < 0.8 at = 2.76 TeV. From these comparisons they were able to 

obtain an estimate of the cold nuclear matter effects on Pb-|-Pb collisions, assuming that the ratios in the forward 
and backward regions factorize. They then used their p+Pb measurements at the higher energy to form the ratio 
'S'j/i/i = RphPh/{RpPh(,+y) X RpPh(,—y))- From this comparison, they were able to deduce that, at least for their 
dimuon measurement, the Pb+Pb data are significantly more surpressed than expected for cold matter effects alone 

0- . . n. 

However, the comparison in Ref. [2[ is made for different rapidity regions due to the asymmetric beams in p-|-Pb 
collisions. In addition, the difference in the nucleon-nucleon center of mass energy is almost a factor of two. These 
two effects will cause the x values probed in p+Pb and Pb+Pb collisions to be shifted so that the correspondence is 
not exact. 

In this section, we compare the RpPh and i?pbPb ratios at LO (rapidity only) and NLO (rapidity andpr in symmetric 
forward and backward regions) at the same energy for p+Pb and Pb+Pb so that the correspondence between the 
X values should be exact. In addition, since we employ the same there no rapidity shift in p+Pb relative to 

Pb+Pb. We can thus check this factorization does indeed explicitly hold for shadowing effects in the GEM. In our 
calculations here, we employ the central EPS09 sets. The results should be similar for other nPDE sets. 
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FIG. 21: (Color online) The mass and scale uncertainties in the T ratio RpPh{pT) in the ALICE acceptance at forward (a), 
backward (b) and central (c) rapidity. The ratio Rpph{y) is compared to the ALICE Q] and LHCb data in (d). The EPS09 
NLO uncertainty band is shown in red while the uncertainties calculated with method vl in blue, v2 in magenta and v3 in 
black. 


It is straightforward to make this comparison at leading order in the color evaporation model (CEM) where the pt 
of the QQ pair is zero and the xi and X 2 values are related to the quarkonium rapidity by xi,X 2 = exp(± 2 /). 

As long as factorization is assumed for the parton densities, the individual shadowing ratios applied to each also 
factorize. The result is compared for the EPS09 LO central shadowing set in Figs. I^ a) and E^ a) for ^/snn = 2.76 
TeV. In these 2 —>• 1 kinematics, the relation RAA(y) = RpA{y) x RpA{—y) is exact. 

Here the dashed and dot-dashed curves show the ratios RpPh at positive and negative rapidities respectively. (In 
this case, the —y refers to the ion beam moving toward positive rapidity while -\-y is the standard assumption that the 
proton beam moves toward positive rapidity, the convention adopted throughout this paper.) The blue dots represent 
the product of i?pPb(+2/) and i?pPb(—2/)- At LO, they lie exactly on top of the red solid curve calculated for Pb-bPb. 

A consequence of the lower energy is that the antishadowing peak moves close to midrapidity. This effect is most 
clearly seen for the T in Fig. [24l a) where the drop into the EMC region and the subsequent rise into the Fermi motion 
region at high x is clearly visible. 
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FIG. 22: (Color online) The mass and scale uncertainties in the T forward-backward ratio Rfb{pt) in the ALICE overlap 
region (a) and Rpaiy) (b). The EPS09 NLO uncertainty band is shown in red while the uncertainties calculated with method 
vl in blue and u3 in black. The ALICE and LHCb data are also shown in (b). 


It is also obvious, both for the J/r/^ and T, that the Pb-|-Pb ratio gives stronger shadowing at midrapidity than at 
forward (and backward) rapidity since, at y ^ 0, the ratios i?ppb at ±y are both less than unity while, at 2.5 < \y\ < 4.0, 
these ratios are in the antishadowing region. Therefore, the combination of the two in Pb-|-Pb interactions is always 
less than unity with stronger shadowing at midrapidity. 

This was also the case for RHIC cold matter A A calculations at = 200 GeV [2^. More suppression is 

predicted at y = 0 than at forward and backward rapidities with all the shadowing parameterizations for both J/tp 
and T production. At RHIC, the A + A data are more suppressed at forward rapidity than at central rapidity, both 
in the minimum bias data as a function of rapidity and as a function of collision centrality, as quantified by the 
number of participant nucleons. Standard models of shadowing alone or shadowing with absorption by nucleons in 
cold nuclear matter or shadowing combined with dissociation in a quark-gluon plasma leads to strong suppression at 
central rapidities. However, J/\p regeneration by dynamical coalescence of c and c quarks in the medium 1^ is 
biased toward central rapidities and could lead to more suppression at forward rapidity relative to central rapidity 
since the rapidity distribution of J/tp production by coalescence is expected to be narrower than the initial J/ip 
rapidity distribution |6ll|. Thus, with coalescence, there should be more suppression at forward y than at midrapidity. 
Statistical coalescence [62|, which is based only on the energy density of the medium, which is higher at y ^ 0, would 
also support this picture. 

The same trend has been observed at the LHC, see Ref. and references therein. Coalescence production of 
the J/ip should be more important than at RHIC since more cc pairs are created in a central Pb-fPb collision at 
^/spppp = 2.76 TeV. We can also expect T production by coalescence at the LHC to be similar to that expected for 
the J/Ip Sit RHIC since the bb production cross section at the LHC will be similar to the cc production cross section 
at RHIC [H. 

At next-to-leading order, the assumption of factorization of shadowing is less straightforward because quarkonium 
production in the CEM includes a large contribution from 2 —?> 3 diagrams so that the correlation between the initial 
momentum fractions xi and X 2 with the rapidity of the quarkonium state is weaker, particularly for high pT- However, 
factorization is seen to still hold as a function of rapidity at NLO, as shown in Figs. E5I b) and IMI b). calculated with 
the EPS09 NLO central set, also at y/sTm = 2.76 TeV. There are some small fluctuations in the points relative to 
the curve but these are within the size of the points and are therefore negligible. 

Note, however, that there is a significant difference in the shape of the shadowing ratios as a function of rapidity 
between the LO and NLO results of Figs. I^ a) and (b) as well as between Figs. I^ a) and (b). As discussed in 
Sec. IVI C[ this is due to the differences in the EPS09 LO and NLO sets themselves and gives some indication of the 
uncertainty inherent in the extraction of the nuclear gluon density. 

Finally, we turn to the pt dependence at NLO, unavailable in the CEM at LO. We might expect to see the largest 
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FIG. 23: (Color online) The J/i/> Raa (red) ratio is compared to the product R-pA{+y) x RpA{—y) (points) along with the 
individual pA ratios at forward (dashed) and backward (dot-dashed) rapidity. Results are compared for the rapidity distributions 
at LO (a) and NLO (b) as well as for the pt dependence at NLO (c). 


deviations in this comparison because of the different kinematics in 2 —?> 2 and 2 —?> 3 interactions in the CEM at 
NLO. At high pp, the 2 —>■ 3 kinematics is dominant. However, the agreement between the factorized product and 
the direct Pb-bPb calculation is very good for both J/'f/; and T up to quite high pp, see Figs. IMT c) and lMT c). (The 
result is shown up to pp = 30 GeV.) There are more fluctuations in the T calculations. Despite this, the agreement 
is still very good. 

The comparisons shown here, at the same energy for both p-|-Pb and Pb-|-Pb collisions, demonstrate that the cold 
matter effects due to shadowing in heavy-ion collisions can be effectively deduced from proton-nucleus collisions, 
preferably at the same energy. 
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FIG. 24: (Color online) The T Raa (red) ratio is compared to the product RpA{+y)xRpA{—y) (points) along with the individual 
pA ratios at forward (dashed) and backward (dot-dashed) rapidity. Results are compared for the rapidity distributions at LO 
(a) and NLO (b) as well as for the pr dependence at NLO (c). 


G. Comparison to RHIC results 

Finally, we compare the EPS09 NLO calculations to the RHIC J/ij) and T data in Figs. [55] and [501 No absorption 
is included here, even though the absorption cross section may be non-negiglble at RHIC energies. The EPS09 
LO calculations were employed by the PHENIX Collaboration to extract the putative absorption cross section 
required to make a shadowing and absorption scenario agree with the data. The NLO calculations are presented as 
a function of pr and rapidity together for the first time. 

The J/ip pT-dependent data, shown in Fig. ISsT al-fc). follows approximately the same trend at forward, backward 
and midrapidity. At pp ~ 0, i?dAu ^ 0.7 — 0.8. The ratio then increases with pp until becoming compatible with unity, 
albeit it with rather poor statistics, at pp > 4 GeV. At forward and midrapidity the calculations agree with the trends 
of the data rather well. At low pp^ the x range is in the shadowing region but moves toward the antishadowing region 
as Pp increases, see the shape of RdKuiy) in Fig. ISST d). However, at backward rapidity at RHIC, the calculations 
suggest that, at low pp, the antishadowing region is probed while, at higher pp, the EMC region is reached. Thus, 
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FIG. 25: (Color online) The J/iji ratio RdAviipr) from PHENIX at forward (a), backward (b) and central (c) rapidity [^. The 
ratio RdAu{y) [HI is shown (d). The EPS09 NLO uncertainty band is compared to the data. 


for Pt < 2 GeV, the trend of the data and the calculations are opposite. For px > 2 GeV, the uncertainties in the 
data and the calculations are large enough for the results to be compatible. 

A recent calculation studied the centrality dependence of gluon shadowing. The impact parameter dependence, 
together with a rapidity-dependent absorption cross section, was extracted from the centrality dependence of RdAuiu) 
The large absorption cross section required at backward rapidity was explained in the context of an expanding 
color octet that reaches its final-state size at backward rapidity but has a negligible effect at forward rapidity because 
the primordial J/ij) passes through the target before reaching its full size |55| . Presumably, this would also have 
the desired effect on the px dependence in the backward region since low px J/ip^s will be strongly absorbed in this 
region while those at higher px will again pass through the transverse direction of the target before fully forming, 
maintaining the agreement of the calculation with the data at higher px- 

Finally, we discuss the T results shown in Fig. [211 There is only limited data for the T, including the low statistics 
PHENIX and STAR data presented in Fig. |26|d). The uncertainties on the measurements, along with the 
large uncertainties in the calculation at backward rapidity, allow the calculation and the data to agree within the 
errors, except for the STAR point at y ^ 0. Although it has the lowest statistical uncertainty, it exhibits strong 
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FIG. 26: (Color online) The T ratio RdAuipr) for RHIC at forward (a), backward (b) and central (c) rapidity. The ratio 
RdAu{y) is also shown (d). The EPS09 NLO uncertainty band is compared to the PHENIX and STAR data on 
RdAuiv) data in (d). 


suppression in a rapidity region where antishadowing is predicted. The px dependence has not yet been measured. 

The calculated px dependence in the forward region, Fig. [26r a). is already entering the antishadowing region for 
px > 6 GeV. At midrapidity, the entire px range is in this region so that RdAuiPr) > 1 for all px- On the other hand, 
at backward rapidity, most of the px region is, in fact, in the EMC region where the uncertainies in the calculated 
RdAniPr) are large, see Fig. I^b). 

It is unlikely that the situation for T can be improved with the current RHIC detectors. However, the future 
sPHENIX detector is expected to be able to separate the three T(S) states. This improved mass resolution, 
together with the higher efficiency expected for sPHENIX and higher luminosity of the sPHENIX runs, means that 
the Px dependence could be measured, if a d+Au run is made during the sPHENIX lifetime. 
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VII. SUMMARY 

The results with EPS09 NLO shadowing agree with the measured RpPh within the uncertainties of both the 
calculation and the data. However, the forward-backward ratios, independent of the interpolated p + p baseline, are 
not as well described. The older EKS98 LO parameterization, although not applied consistently in a NLO calculation, 
does the best job of describing all the data. 

The shadowing parameterizations used in our study exhibit a wide range of behavior for the nuclear gluon density 
at low X, an x region outside the current range of the fits from fixed-target nDIS data at higher x and low If 
nuclear data were available from high energy e + A collisions, the nuclear gluon densities could be more precisely 
pinned down by global analyses of the scale dependence of the nuclear structure functions. In hadroproduction, direct 
photon or open charm production, dominated by gluon-induced processes but without the additional complexities 
of nuclear absorption, could be utilized to study the nuclear gluon density. The LHC data, with the lower x reach, 
should be incorporated into a new fit of the nuclear parton densities. 

We note that, since we have assumed absorption is negligible at the LHC and include no other cold nuclear matter 
effect, the uncertainties on the ratios can be obtained from the EPS09 NLO bands shown in the figures. However, if 
other effects are incorporated, a more extensive error analysis, including the uncertainties on these other effects, is 
necessary. 
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